function generateMfuncResidualExpectedControl_file(allModelParams,x,y,xp,yp,nx1,ne,h11,heteroShock,positionFile,nameOfFunction)
% This function generates an M function which is the objective function
% when computing the expected value of any control variable by projection

% We start deleting the old version of the file - if it exists
saveFile = [positionFile,'\',nameOfFunction];
if exist(saveFile,'file') > 0
    delete(saveFile)
end

text = ['function [res,resEuler,Evar_cu] = ',nameOfFunction(1:end-2),'(thetaVec,setup)'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% This function computes the residuals for the problem for computing E_t[var_cup] ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% We use the notation that Evar_cu denotes E_t[var_cup], meaning that the';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% problem is -Evar_cu + var_cup = 0';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Getting some information from setup ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'scaleX_in_g = setup.scaleX_in_g;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'out         = setup.out;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'x_cu        = setup.xGrid;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'numX        = size(x_cu,1);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nNodes      = setup.n_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'wNodes      = transpose(setup.weight_nodes);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'epsNodes    = setup.epsi_nodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'params      = setup.params;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'appOrder    = setup.appOrder;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx          = setup.nx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'ny          = setup.ny;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'nx1         = setup.nx1;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gSS         = setup.gSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'hSS         = setup.hSS;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The unknowns: for Evar_cu = E_t[var_cup]';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'thetaProj   = reshape(thetaVec,1,setup.nDim2Theta);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'g0          = thetaProj(1,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'gx          = thetaProj(1,2:1+nx); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxV    = thetaProj(1,nx+2:1+nx+nx*(nx+1)/2); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    gxxxV   = thetaProj(1,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The knowns: for var_cu  ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'thetaVar    = reshape(setup.thetaVar,1,setup.nDim2Theta);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'var0        = thetaVar(1,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'varx        = thetaVar(1,2:1+nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    varxx   = thetaVar(1,nx+2:1+nx+nx*(nx+1)/2);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2 ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '    varxxx  = thetaVar(1,1+nx+nx*(nx+1)/2+1:1+nx+nx*(nx+1)/2+nx*(nx+1)*(nx+2)/6); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Unfold params ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(allModelParams)
    text = [allModelParams{i}, '= params.',allModelParams{i},';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end 
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% Unfold the states and controls for the core of the model - all in levels ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the states at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
     text = [char(x(i)),' = out.',char(x(i)),';'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(y)
    text = [char(y(i)),' = out.',char(y(i)),';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the states at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(xp)
     text = [char(xp(i)),' = out.',char(xp(i)),';'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Getting the controls at time t+1 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(yp)
     text = [char(yp(i)),' = out.',char(yp(i)),';'];
     dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Computing the new control variable  ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    text = ['x',num2str(i),'_cu ',' = x_cu(:,',num2str(i),')/scaleX_in_g(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'Egfunc_cu = g0(1,1)';
for j=1:length(x)
    name = [name,['+gx(1',',',num2str(j),').*x',num2str(j),'_cu']];
end
text = [name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
idx = 0;
name = 'Egfunc_cu = Egfunc_cu';
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        idx = idx + 1;
        if alfa1 == alfa2
            name = [name,['+gxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
        else
            name = [name,['+2.*gxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu']];
        end
    end
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'Egfunc_cu = Egfunc_cu';
idx = 0;
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        for alfa3=alfa2:length(x)
            idx = idx + 1;
            if alfa1 == alfa2 && alfa2 == alfa3
                name = [name,['+gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            elseif alfa1 == alfa2 && alfa2 ~= alfa3
                name = [name,['+3.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            elseif alfa1 ~= alfa2 && alfa2 == alfa3
                name = [name,['+3.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                name = [name,['+6.*gxxxV(1',',',num2str(idx),').*x',num2str(alfa1),'_cu','.*x',num2str(alfa2),'_cu','.*x',num2str(alfa3),'_cu']];
            else
                [alfa1,alfa2,alfa3];
            end
        end
    end
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%% The control variable in the next period ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xVec_cup = zeros(numX,nx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'mx = setup.mx;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if mx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
nx11 = length(h11);
for i=1:nx11
    TEMP = char(h11(i,:));
    TEMP = strrep(TEMP,'*','.*');
    TEMP = strrep(TEMP,'/','./');
    TEMP = strrep(TEMP,'^','.^');
    text = ['    ',char(xp(1,i)),' = ',TEMP,';'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if setup.myx > 0 ';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1-nx11
    text = ['    xVec_cup(:,mx+',num2str(i),') = ',char(y(1,i)),'(:,1)-gSS(',num2str(i),',1);'];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = 'end ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'diaghx = diag(setup.hx);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'xVec_cup(:,nx1+1:nx)  = x_cu(:,nx1+1:nx).*repmat(transpose(diaghx(nx1+1:nx)),numX,1); ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '  ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '%Adding shocks to the states at t + 1';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'diagEta = diag(setup.eta(nx1+1:nx,:));';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:nx1
    name = ['x',num2str(i),'_cup  = repmat(xVec_cup(:,',num2str(i),'),1,nNodes);'];
    text = name;
    dlmwrite(saveFile,text,'-append','delimiter','');
end
for i=1:ne
    name = ['x',num2str(nx1+i),'_cup  = repmat(xVec_cup(:,',num2str(nx1+i),'),1,nNodes) +',...
         '2./(exp(-',char(heteroShock(i,1)),'.*x_cu(:,',num2str(nx1+i),'))+1).*diagEta(',...
         num2str(i),',1).*repmat(epsNodes(',num2str(i),',:),numX,1);'];
    text = ['',name];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Scaling the states at time t+1 by the scaleX_in_g';
dlmwrite(saveFile,text,'-append','delimiter','');
for i=1:length(x)
    name = ['x',num2str(i),'_cup  = x',num2str(i),'_cup/scaleX_in_g(',num2str(i),',1);'];
    text = ['',name];
    dlmwrite(saveFile,text,'-append','delimiter','');
end
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'var_cup = var0(1,1)';
for j=1:length(x)
    name = [name,['+varx(1',',',num2str(j),').*x',num2str(j),'_cup']];
end
text = [name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 1';
dlmwrite(saveFile,text,'-append','delimiter','');
name = 'var_cup = var_cup';
idx = 0;
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        idx = idx + 1;
        if alfa1 == alfa2
            name = [name,['+varxx(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
        else
            name = [name,['+2.*varxx(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup']];
        end
    end
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'if appOrder > 2';
dlmwrite(saveFile,text,'-append','delimiter','');
idx = 0;
name = 'var_cup = var_cup';
for alfa1=1:length(x)
    for alfa2=alfa1:length(x)
        for alfa3=alfa2:length(x)
            idx = idx + 1;
            if alfa1 == alfa2 && alfa2 == alfa3
                name = [name,['+varxxx(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            elseif alfa1 == alfa2 && alfa2 ~= alfa3
                name = [name,['+3.*varxxx(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            elseif alfa1 ~= alfa2 && alfa2 == alfa3
                name = [name,['+3.*varxxx(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            elseif alfa1 ~= alfa2 && alfa2 ~= alfa3
                name = [name,['+6.*varxxx(1',',',num2str(idx),').*x',num2str(alfa1),'_cup','.*x',num2str(alfa2),'_cup','.*x',num2str(alfa3),'_cup']];
            else
                [alfa1,alfa2,alfa3];
            end
        end
    end
end
text = ['    ',name,';'];
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
text = ' ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% Numerical integration';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'Evar_cu = var_cup*wNodes;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = '% The Euler errors ';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'resEuler = -Egfunc_cu + Evar_cu;';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'resEulerWeight = resEuler.*(1./(1+sum(x_cu.^2,2)));';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'res = resEulerWeight(:);';
dlmwrite(saveFile,text,'-append','delimiter','');
text = 'end';
dlmwrite(saveFile,text,'-append','delimiter','');
end